load packages and functions
source('mytheme.R')
# model version
modelversion = 'gap_log_rstan'
rstanmodelPath = 'modelrlt'
modelPath = paste0(rstanmodelPath, '/models/', modelversion)
Merge the model Result data
load data
AllExpData = read.csv(paste0("../data/AllValidData.csv"))
dur <- sort(unique(AllExpData$curDur))
AllExpData$WMSize <- factor(AllExpData$WMSize, labels = c("low", "medium", "high"))
# 1: 500ms, 2: 2500, 3 : 2000ms the mean reaction time of WM test
AllExpData$gap <- factor(AllExpData$gap, labels = c('short','long', 'not sure'))
AllExpData[which(AllExpData$Exp == 'Exp4'),"Exp"] = "Exp4a"
AllExpData[which(AllExpData$Exp == 'Exp5'),"Exp"] = "Exp4b"
Corrct rate
WMCrr2 <- ggplot(meanForPlot, aes(WMSize, mean_WMCrr, ymin = mean_WMCrr - se_WMCrr, ymax = mean_WMCrr + se_WMCrr,
group =Exp, color = Exp, fill = Exp))+
geom_line(stat = "identity",position = position_dodge(width = 0.2))+
geom_point(stat = "identity",position = position_dodge(width = 0.2))+
geom_errorbar(width=.2, position = position_dodge(width = 0.2)) +
coord_cartesian(ylim = c(0.5, 1)) +
colorSet5+
labs(x = "Memory load", y = "Mean accuracy in WM task") +
theme_new
WMCrr2
ggsave(paste0(getwd(), "/figures/WMCrr2.png"), WMCrr2, width = 4, height = 4)
### generate WM correct rates
AllExpData$WMCrr <- AllExpData$TPresent == AllExpData$WMRP
m_wmp<- dplyr::group_by(AllExpData, Exp, WMSize, NSub) %>%
dplyr::summarize(m_WMCrr = mean(WMCrr), n =n(), se_WMCrr = sd(WMCrr)/sqrt(n-1))
## `summarise()` has grouped output by 'Exp', 'WMSize'. You can override using the
## `.groups` argument.
comparison between Exp.4a and Exp.4b)
correct rate
#plot WM correct rates
dplyr::group_by(AllExpData%>%filter(Exp %in% c('Exp4a', 'Exp4b')), Exp, WMSize, NSub, gap) %>%
dplyr::summarize(m_WMCrr = mean(WMCrr), n = n(), se_WMCrr = sd(WMCrr)/sqrt(n-1)) -> correctrate_gap
## `summarise()` has grouped output by 'Exp', 'WMSize', 'NSub'. You can override
## using the `.groups` argument.
write.csv(correctrate_gap, paste0(modelPath, '/rlt/correctrate_gap.csv'))
correctrate_gap%>%
dplyr::group_by(Exp, WMSize, gap)%>%
dplyr::summarize( n = n(),
mean_WMCrr = mean(m_WMCrr), se_WMCrr = sd(m_WMCrr)/sqrt(n-1) ) -> meanForPlot_gap
## `summarise()` has grouped output by 'Exp', 'WMSize'. You can override using the
## `.groups` argument.
WMCrr_gap <- ggplot(meanForPlot_gap, aes(WMSize, mean_WMCrr, ymin = mean_WMCrr - se_WMCrr, ymax = mean_WMCrr + se_WMCrr, group =interaction(Exp,gap), color = gap, shape = Exp, linetype = Exp))+
geom_line(stat = "identity",position = position_dodge(width = 0.2))+
geom_point(stat = "identity",position = position_dodge(width = 0.2))+
geom_errorbar(width=.2, position = position_dodge(width = 0.2)) +
coord_cartesian(ylim = c(0.5, 1)) +
colorSet5+
labs(x = "Memory load", y = "Mean accuracy in WM task") +theme_new
WMCrr_gap
ggsave(paste0(getwd(), "/figures/WMCrr_gap.png"), WMCrr_gap, width = 4, height = 4)
observed RP
AllExpData%>% filter(Exp %in% c('Exp4a', 'Exp4b')) %>%
dplyr::group_by(Exp, curDur, WMSize, NSub, gap) %>%
dplyr::summarize(n = n(),
m_repDur = mean(repDur),
se_repDur = sd(repDur)/sqrt(n-1)) ->RP_obs_gap_bias
## `summarise()` has grouped output by 'Exp', 'curDur', 'WMSize', 'NSub'. You can
## override using the `.groups` argument.
RP_obs_gap_bias$rep_bias = RP_obs_gap_bias$m_repDur - RP_obs_gap_bias$curDur
ezANOVA(data = RP_obs_gap_bias%>%filter(Exp =='Exp4b'), dv= rep_bias, wid=NSub, within= .(curDur, gap, WMSize) )
## Warning: Converting "NSub" to factor for ANOVA.
## Warning: "curDur" will be treated as numeric.
## Warning: You have removed one or more levels from variable "gap". Refactoring
## for ANOVA.
## Warning: There is at least one numeric within variable, therefore aov() will be
## used for computation and no assumption checks will be obtained.
## $ANOVA
## Effect DFn DFd F p p<.05 ges
## 1 curDur 1 15 96.33440029 6.384821e-08 * 0.8176127533
## 2 gap 1 15 0.03549505 8.530873e-01 0.0001018033
## 3 WMSize 2 30 9.48732989 6.415750e-04 * 0.0895670700
## 4 curDur:gap 1 15 16.35270907 1.060464e-03 * 0.0220610259
## 5 curDur:WMSize 2 30 9.34453850 7.003993e-04 * 0.0280287517
## 6 gap:WMSize 2 30 2.51037940 9.816022e-02 0.0027254302
## 7 curDur:gap:WMSize 2 30 1.63535038 2.117823e-01 0.0021874657
RP_jasp <- RP_obs_gap_bias%>%filter(Exp =='Exp4b')%>% select(c( "curDur", "WMSize","NSub","gap", "rep_bias"))%>%
pivot_wider(names_from = c("curDur", "WMSize", "gap"), values_from = rep_bias, names_sep="")
## Adding missing grouping variables: `Exp`
#write.csv(RP_jasp, 'my_bias_rep_Exp5.csv')
RP_obs_gap <- ggplot( data = RP_obs_gap_bias%>%
dplyr::group_by(Exp, curDur, WMSize, gap) %>%
dplyr::summarize(m_m_repDur = mean(m_repDur),
m_se_repDur = mean(se_repDur)), aes(x = curDur, y = m_m_repDur, color=as.factor(WMSize), shape = gap)) +
geom_point()+
geom_line()+
geom_abline(slope=1, intercept=0)+
facet_grid(cols = vars(Exp)) +
labs(x="Sample intervals (s)", y="Reproduction (s)", shape=" ", color = "Memory Load")+
theme_new+colorSet3+
scale_x_continuous(breaks=seq(0, 1.6, 0.4))+ theme(legend.position="top")
## `summarise()` has grouped output by 'Exp', 'curDur', 'WMSize'. You can override
## using the `.groups` argument.
ggsave(paste0(getwd(), "/", modelPath, "/figures/RP_obs.png"), RP_obs_gap, width = 6, height = 4)
RP_obs_gap
central tendency effect (slope) and IPs (linear regression)
#Observed Indifference Point
obs_model <- function(df) {
lm(repDur ~ curDur, data = df)
}
obs_InP <- AllExpData %>%
dplyr::group_by(NSub, Exp, WMSize, gap) %>% nest() %>%
mutate(model = map(data, obs_model)) %>% # linear regression
mutate(slope = map(model, broom::tidy)) %>% # get estimates
unnest(slope, .drop = TRUE) %>% # remove raw data
select(-std.error,-statistic, -p.value) %>% # remove unnessary columns
spread(term, estimate) %>% # spread stimates
dplyr::rename(Intercept = `(Intercept)`, slope = curDur) # rename columns
## Warning: The `.drop` argument of `unnest()` is deprecated as of tidyr 1.0.0.
## All list-columns are now preserved.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was generated.
obs_InP$model = NULL
obs_InP$data = NULL
obs_InP$inP = obs_InP$Intercept /(1-obs_InP$slope)
obs_InP_Exp1_4 <- AllExpData %>%
dplyr::group_by(NSub, Exp, WMSize) %>% nest() %>%
mutate(model = map(data, obs_model)) %>% # linear regression
mutate(slope = map(model, broom::tidy)) %>% # get estimates
unnest(slope, .drop = TRUE) %>% # remove raw data
select(-std.error,-statistic, -p.value) %>%
spread(term, estimate) %>% # spread stimates
dplyr::rename(Intercept = `(Intercept)`, slope = curDur) # rename columns
obs_InP_Exp1_4$model = NULL
obs_InP_Exp1_4$data = NULL
obs_InP_Exp1_4$inP = obs_InP_Exp1_4$Intercept /(1-obs_InP_Exp1_4$slope)
obs_InP_slope_Exp1_jasp <- obs_InP_Exp1_4 %>% filter(Exp == 'Exp1') %>% select(c("WMSize", "NSub","inP", "slope"))%>% pivot_wider(names_from = c("WMSize"), values_from = c(inP, slope), names_sep="_")
## Adding missing grouping variables: `Exp`
write.csv(obs_InP_slope_Exp1_jasp, paste0(modelPath, '/rlt/obs_InP_slope_Exp1_jasp.csv'))
obs_InP_slope_Exp2_jasp <- obs_InP_Exp1_4 %>% filter(Exp == 'Exp2') %>% select(c("WMSize", "NSub","inP", "slope"))%>% pivot_wider(names_from = c("WMSize"), values_from = c(inP, slope), names_sep="_")
## Adding missing grouping variables: `Exp`
write.csv(obs_InP_slope_Exp2_jasp, paste0(modelPath, '/rlt/obs_InP_slope_Exp2_jasp.csv'))
obs_InP_slope_Exp3_jasp <- obs_InP_Exp1_4 %>% filter(Exp == 'Exp3') %>% select(c("WMSize", "NSub","inP", "slope"))%>% pivot_wider(names_from = c("WMSize"), values_from = c(inP, slope), names_sep="_")
## Adding missing grouping variables: `Exp`
write.csv(obs_InP_slope_Exp3_jasp, paste0(modelPath, '/rlt/obs_InP_slope_Exp3_jasp.csv'))
obs_InP_slope_Exp4a_jasp <- obs_InP_Exp1_4 %>% filter(Exp =='Exp4a') %>% select(c("WMSize", "NSub", "inP","slope"))%>%
pivot_wider(names_from = c("WMSize"), values_from = c(inP, slope), names_sep="_")
## Adding missing grouping variables: `Exp`
write.csv(obs_InP_slope_Exp4a_jasp, paste0(modelPath, '/rlt/obs_InP_slope_Exp4a_jasp.csv'))
obs_InP_Exp4b_jasp <- obs_InP %>% filter(Exp == 'Exp4b') %>% select(c("WMSize", "NSub", "gap","inP", "slope"))%>% pivot_wider(names_from = c("WMSize", "gap"), values_from = c(inP,slope), names_sep="_")
## Adding missing grouping variables: `Exp`
write.csv(obs_InP_Exp4b_jasp, paste0(modelPath, '/rlt/obs_InP_Exp4b_jasp.csv'))
mobs_InP = obs_InP %>% group_by(Exp, WMSize, gap)%>%
dplyr::summarise(n=n(),
m_Intercept = mean(Intercept),
se_Intercept= sd(Intercept)/sqrt(n-1),
m_inP = mean(inP),
se_inP = sd(inP)/sqrt(n-1),
m_slope = mean(slope),
se_slope = sd(slope)/sqrt(n-1))
## `summarise()` has grouped output by 'Exp', 'WMSize'. You can override using the
## `.groups` argument.
plt_InP_linear_gap<- ggplot(data = mobs_InP, aes(x=WMSize, y=m_inP, group = gap, color = gap))+
geom_line(stat = "identity",position = position_dodge(width = 0.2))+
geom_point(stat = "identity",position = position_dodge(width = 0.2))+
geom_errorbar(width=.3, aes(ymin = m_inP - se_inP, ymax = m_inP + se_inP), position = position_dodge(width = 0.2)) +theme_new+
labs(colour = "Gap")+colorSet3+
facet_wrap(~Exp)+
xlab(' ')+ylab("indifference point (s)")+guides(shape="none")+
theme(legend.position = "top")
ggsave(paste0(getwd(), "/", modelPath, "/figures/plt_InP_linear_gap.png"), plt_InP_linear_gap, width = 5, height = 5)
plt_InP_linear_gap
ezANOVA(data = obs_InP%>%filter(Exp =='Exp4b'), dv= inP, wid=NSub, within= .(gap, WMSize) )
## Warning: Converting "NSub" to factor for ANOVA.
## Warning: You have removed one or more levels from variable "gap". Refactoring
## for ANOVA.
## $ANOVA
## Effect DFn DFd F p p<.05 ges
## 2 gap 1 15 4.065468 0.062041074 0.019543220
## 3 WMSize 2 30 8.762649 0.001006806 * 0.078253848
## 4 gap:WMSize 2 30 0.584382 0.563670241 0.003061788
##
## $`Mauchly's Test for Sphericity`
## Effect W p p<.05
## 3 WMSize 0.4515805 0.003829536 *
## 4 gap:WMSize 0.5210587 0.010428128 *
##
## $`Sphericity Corrections`
## Effect GGe p[GG] p[GG]<.05 HFe p[HF] p[HF]<.05
## 3 WMSize 0.6458198 0.004991546 * 0.6808339 0.004255437 *
## 4 gap:WMSize 0.6761593 0.502620038 0.7194299 0.512209537
plt_RP_slope_linear_gap<- ggplot(data = mobs_InP, aes(x= WMSize, y=m_slope, group = gap,color = gap, shape = gap))+
geom_line(stat = "identity",position = position_dodge(width = 0.2))+
geom_point(stat = "identity",position = position_dodge(width = 0.2))+
geom_errorbar(width=.3, aes(ymin = m_slope - se_slope, ymax = m_slope + se_slope), position = position_dodge(width = 0.2)) +theme_new+
facet_wrap(~Exp)+
labs(colour = "Gap", shape = "Gap")+colorSet3+
xlab(' ')+ylab("Slope of reproduction (s)")+
theme(legend.position = "top")
ggsave(paste0(getwd(), "/", modelPath, "/figures/plt_RP_slope_linear_gap.png"), plt_RP_slope_linear_gap, width = 5, height = 5)
plt_RP_slope_linear_gap
ezANOVA(data = obs_InP%>%filter(Exp =='Exp4b'), dv= slope, wid=NSub, within= .(gap, WMSize) )
## Warning: Converting "NSub" to factor for ANOVA.
## Warning: You have removed one or more levels from variable "gap". Refactoring
## for ANOVA.
## $ANOVA
## Effect DFn DFd F p p<.05 ges
## 2 gap 1 15 16.352709 0.0010604639 * 0.027930881
## 3 WMSize 2 30 9.344539 0.0007003993 * 0.035428935
## 4 gap:WMSize 2 30 1.635350 0.2117822688 0.002784549
##
## $`Mauchly's Test for Sphericity`
## Effect W p p<.05
## 3 WMSize 0.8669475 0.36808622
## 4 gap:WMSize 0.6136631 0.03277255 *
##
## $`Sphericity Corrections`
## Effect GGe p[GG] p[GG]<.05 HFe p[HF] p[HF]<.05
## 3 WMSize 0.8825716 0.001234014 * 0.9914082 0.0007299619 *
## 4 gap:WMSize 0.7213254 0.219115993 0.7775271 0.2179558985
load model results
check model parameters
### Average Parameters
mm_Baypar <- dplyr::group_by(Bayparlist, Exp) %>%
dplyr::summarize( m_sig_s2 = mean(sig_s2),
m_sig_pr2_log = mean(sig_pr2_log),
m_ks= mean(ks),
m_kr = mean(kr),
m_ls = mean(ls),
m_ts = mean(ts),
m_mu_pr = mean(mu_pr),
m_sig_pr2 = mean(sig_pr2),
m_mu_pr_log= mean(mu_pr_log),
m_sig_mn2 = mean(sig_mn2),
n= n(),
se_sig_s2 = sd(sig_s2)/sqrt(n-1),
se_sig_mn2 = sd(sig_mn2)/sqrt(n-1),
se_sig_pr2_log = sd(sig_pr2_log)/sqrt(n-1),
se_ks = sd(ks)/sqrt(n-1),
se_kr = sd(kr)/sqrt(n-1),
se_ls =sd(ls)/sqrt(n-1),
se_mu_pr_log = sd(mu_pr_log)/sqrt(n-1))
mm_Baypar
Average Parameters
mBaypar_sub <- dplyr::group_by(Bayparlist, Exp, NSub) %>%
dplyr::summarize( m_sig_s2 = mean(sig_s2),
m_sig_pr2_log = mean(sig_pr2_log),
m_ks= mean(ks),
m_kr = mean(kr),
m_ls = mean(ls),
m_ts = mean(ts),
m_mu_pr = mean(mu_pr),
m_sig_pr2 = mean(sig_pr2),
m_mu_pr_log= mean(mu_pr_log),
n= n(),
se_sig_s2 = sd(sig_s2)/sqrt(n-1),
se_sig_mn2 = sd(sig_mn2)/sqrt(n-1),
se_sig_pr2_log = sd(sig_pr2_log)/sqrt(n-1),
se_ks = sd(ks)/sqrt(n-1),
se_kr = sd(kr)/sqrt(n-1),
se_ls =sd(ls)/sqrt(n-1),
se_mu_pr_log = sd(mu_pr_log)/sqrt(n-1))
## `summarise()` has grouped output by 'Exp'. You can override using the `.groups`
## argument.
fig_ks_subj = ggplot(mBaypar_sub, aes(Exp, m_ks, color = Exp)) +
geom_boxplot(position = position_dodge()) +
geom_jitter(shape=16, position=position_jitter(0.2))+
#facet_wrap(~Exp)+
theme_new+ colorSet5 +
theme(strip.background = element_blank()) + # remove subtitle background
labs(x = ' ', y = 'scale factor Ks')+
theme(legend.position='bottom')
ggsave(paste0(getwd(), "/", modelPath, "/figures/fig_ks_subj.png"), fig_ks_subj, width = 6, height = 6)
fig_ks_subj
fig_kr_subj = ggplot(mBaypar_sub, aes(Exp, m_kr, color = Exp)) +
geom_boxplot(position = position_dodge()) +
geom_jitter(shape=16, position=position_jitter(0.2))+
#facet_wrap(~Exp)+
theme_new+ colorSet5 +
theme(strip.background = element_blank()) + # remove subtitle background
labs(x = ' ', y = 'scale factor Kr')+
theme(legend.position='bottom')
ggsave(paste0(getwd(), "/", modelPath, "/figures/fig_kr_subj.png"), fig_kr_subj, width = 6, height = 6)
fig_kr_subj
fig_ts_subj = ggplot(mBaypar_sub, aes(Exp, m_ts, color = Exp)) +
geom_boxplot(position = position_dodge()) +
geom_jitter(shape=16, position=position_jitter(0.2))+
theme_new+ colorSet5 +
theme(strip.background = element_blank()) + # remove subtitle background
labs(x = ' ', y = 'scale factor ts')+
theme(legend.position='bottom')
ggsave(paste0(getwd(), "/", modelPath, "/figures/fig_ts_subj.png"), fig_ts_subj, width = 6, height = 6)
fig_ts_subj
fig_ls_subj = ggplot(mBaypar_sub, aes(Exp, m_ls, color = Exp)) +
geom_boxplot(position = position_dodge()) +
geom_jitter(shape=16, position=position_jitter(0.2))+
#facet_wrap(~Exp)+
theme_new+ colorSet5 +
theme(strip.background = element_blank()) + # remove subtitle background
labs(x = ' ', y = 'scale factor ls')+
theme(legend.position='bottom')
ggsave(paste0(getwd(), "/", modelPath, "/figures/mfig_ls_subj.png"), fig_ls_subj, width = 6, height = 6)
fig_ls_subj
fig_mu_pr_subj = ggplot(mBaypar_sub, aes(Exp, m_mu_pr, color = Exp)) +
geom_boxplot(position = position_dodge()) +
geom_jitter(shape=16, position=position_jitter(0.2))+
theme_new+ colorSet5 +
theme(strip.background = element_blank()) + # remove subtitle background
labs(x = ' ', y = 'mean of prior')+
theme(legend.position='bottom')
ggsave(paste0(getwd(), "/", modelPath, "/figures/fig_mu_pr_subj.png"), fig_mu_pr_subj, width = 6, height = 6)
fig_mu_pr_subj
fig_sig_pr2_subj = ggplot(mBaypar_sub, aes(Exp, m_sig_pr2, color = Exp)) +
geom_boxplot(position = position_dodge()) +
geom_jitter(shape=16, position=position_jitter(0.2))+
theme_new+ colorSet5 +
theme(strip.background = element_blank()) + # remove subtitle background
labs(x = ' ', y = 'mean of variance')+
theme(legend.position='bottom')
ggsave(paste0(getwd(), "/", modelPath, "/figures/fig_sig_pr2_subj.png"), fig_sig_pr2_subj, width = 6, height = 6)
fig_sig_pr2_subj
Prediction results
#load prediction
AllDat_predY <- read.csv(paste0(modelPath, "/rlt/AllDat_predY_",modelversion,".csv"))
AllDat_predY$WMSize <- as.factor(AllDat_predY$WMSize)
levels(AllDat_predY$WMSize) = c("low", "medium", "high")
AllDat_predY$pred_Bias = AllDat_predY$mu_r - AllDat_predY$curDur
AllDat_predY$predErr = AllDat_predY$mu_r - AllDat_predY$repDur
AllDat_predY$relatErr = AllDat_predY$predErr / AllDat_predY$repDur
AllDat_predY[which(AllDat_predY$Exp == "Exp4"),"Exp"] = "Exp4a"
AllDat_predY[which(AllDat_predY$Exp == "Exp5"),"Exp"] = "Exp4b"
#calculate the mean reproduction biases for the five given intervals for all subjects
mpredY_sub <- dplyr::group_by(AllDat_predY, curDur, Exp, NSub, WMSize) %>%
dplyr::summarize(n = n(),
m_repDur = mean(repDur),
sd_repDur = sd(repDur),
m_mu_r = mean(mu_r),
m_sig_r = mean(sig_r),
m_wp = mean(wp),
se_wp = sd(wp)/sqrt(n-1),
log_lik =mean(log_lik),
cv =sd_repDur/ m_repDur,
pred_cv = mean(sig_r/mu_r),
predRP_err = mean(m_mu_r-m_repDur),
predVar_err = mean(m_sig_r-sd_repDur),
predRP_rerr = mean(abs(m_mu_r-m_repDur)/m_repDur),
predVar_rerr = mean(abs(m_sig_r-sd_repDur)/sd_repDur),
predcv_err = pred_cv-cv,
predcv_rerr = mean(abs(pred_cv-cv)/cv))
## `summarise()` has grouped output by 'curDur', 'Exp', 'NSub'. You can override
## using the `.groups` argument.
mpredY_sub_split <-split(mpredY_sub %>%select(c('Exp', 'NSub', 'WMSize', 'm_repDur', 'curDur')), mpredY_sub$curDur)
mpredY_sub_cv_split <-split(mpredY_sub %>%select(c('Exp', 'NSub', 'WMSize', 'cv', 'curDur')), mpredY_sub$curDur)
mpredY_sub_split1 <-split(mpredY_sub %>%select(c('Exp', 'NSub', 'WMSize', 'curDur' ,'m_repDur')), mpredY_sub$WMSize)
mpredY_sub_cv_split1 <-split(mpredY_sub %>%select(c('Exp', 'NSub', 'WMSize','curDur', 'cv')), mpredY_sub$WMSize)
mpredY_sub_jasp_RP = NULL
mpredY_sub_jasp_cv = NULL
mpredY_sub_jasp_RP1 = NULL
mpredY_sub_jasp_cv1 = NULL
for(i in 1: length(mpredY_sub_split)){
temp = mpredY_sub_split[[i]]
curDur = unique(temp$curDur)
temp$curDur = NULL
colnames(temp) = c('Exp', 'NSub', 'WMSize', paste0('m_repDur_', curDur))
temp_cv = mpredY_sub_cv_split[[i]]
curDur = unique(temp_cv$curDur)
temp_cv$curDur = NULL
colnames(temp_cv) = c('Exp', 'NSub', 'WMSize', paste0('cv_', curDur))
if(i == 1){
mpredY_sub_jasp_RP = temp
mpredY_sub_jasp_cv = temp_cv
}
else{
mpredY_sub_jasp_RP = left_join(mpredY_sub_jasp_RP, temp, by=c("Exp", "NSub", "WMSize"))
mpredY_sub_jasp_cv = left_join(mpredY_sub_jasp_cv, temp_cv, by=c("Exp", "NSub", "WMSize"))
}
}
for (i in 1: length(mpredY_sub_split1)){
temp1 = mpredY_sub_split1[[i]]
WMSize = unique(temp1$WMSize)
temp1$WMSize = NULL
colnames(temp1) = c('Exp', 'NSub', 'curDur', paste0('m_repDur_', WMSize))
temp_cv1 = mpredY_sub_cv_split1[[i]]
WMSize = unique(temp_cv1$WMSize)
temp_cv1$WMSize = NULL
colnames(temp_cv1) = c('Exp', 'NSub', 'curDur', paste0('cv_', WMSize))
if(i == 1){
mpredY_sub_jasp_RP1 = temp1
mpredY_sub_jasp_cv1 = temp_cv1
}
else{
mpredY_sub_jasp_RP1 = left_join(mpredY_sub_jasp_RP1, temp1, by=c("Exp", "NSub", "curDur"))
mpredY_sub_jasp_cv1 = left_join(mpredY_sub_jasp_cv1, temp_cv1, by=c("Exp", "NSub", "curDur"))
}
}
write_csv(mpredY_sub_jasp_RP, paste0(modelPath, '/rlt/RP_Bias_jasp.csv'))
write_csv(mpredY_sub_jasp_cv, paste0(modelPath, '/rlt/cv_jasp.csv'))
write_csv(mpredY_sub_jasp_RP1, paste0(modelPath, '/rlt/RP_Bias_jasp1.csv'))
write_csv(mpredY_sub_jasp_cv1, paste0(modelPath, '/rlt/cv_jasp1.csv'))
write_csv(dplyr::group_by(mpredY_sub, curDur, NSub) %>%
dplyr::summarize(m_cv = mean(cv))%>%spread(curDur, m_cv), paste0(modelPath, '/rlt/m_cv.csv'))
## `summarise()` has grouped output by 'curDur'. You can override using the
## `.groups` argument.
mpredY_sub$RP_bias = mpredY_sub$m_repDur -mpredY_sub$curDur
mpredY_sub_new <- dplyr::group_by(mpredY_sub, curDur, Exp, NSub) %>%
dplyr::summarize(m_RP_bias = mean(RP_bias))%>% spread(curDur, m_RP_bias)
## `summarise()` has grouped output by 'curDur', 'Exp'. You can override using the
## `.groups` argument.
write_csv(mpredY_sub_new%>%filter(Exp == 'Exp1'), paste0(modelPath, '/rlt/RP_Bias_exp1.csv'))
write_csv(mpredY_sub_new%>%filter(Exp == 'Exp2'), paste0(modelPath, '/rlt/RP_Bias_exp2.csv'))
write_csv(mpredY_sub_new%>%filter(Exp == 'Exp3'), paste0(modelPath, '/rlt/RP_Bias_exp3.csv'))
write_csv(mpredY_sub_new%>%filter(Exp == 'Exp4a'), paste0(modelPath, '/rlt/RP_Bias_exp4a.csv'))
write_csv(mpredY_sub_new%>%filter(Exp == 'Exp4b'), paste0(modelPath, '/rlt/RP_Bias_exp4b.csv'))
mpredY_sub_WMsize <- dplyr::group_by(mpredY_sub, WMSize, Exp, NSub) %>%
dplyr::summarize(m_RP_bias = mean(RP_bias))%>% spread(WMSize, m_RP_bias)
## `summarise()` has grouped output by 'WMSize', 'Exp'. You can override using the
## `.groups` argument.
write_csv(mpredY_sub_WMsize%>%filter(Exp == 'Exp3'), paste0(modelPath, '/rlt/RP_Bias_WMsize_exp3.csv'))
write_csv(mpredY_sub_WMsize%>%filter(Exp == 'Exp4a'), paste0(modelPath, '/rlt/RP_Bias_WMsize_exp4a.csv'))
write_csv(mpredY_sub_WMsize%>%filter(Exp == 'Exp4b'), paste0(modelPath, '/rlt/RP_Bias_WMsize_exp4b.csv'))
#### predicted data
m_predY <- mpredY_sub%>%
dplyr::group_by(Exp, curDur, WMSize) %>%
dplyr::summarize(m_m_repDur = mean(m_repDur),
m_sd_repDur = mean(sd_repDur),
m_m_sig_r =mean(m_sig_r),
m_m_mu_r = mean(m_mu_r),
m_m_wp = mean(m_wp),
n = n(),
m_se_wp = sd(se_wp)/sqrt(n-1),
log_lik =mean(log_lik),
mpredRP_err = mean(predRP_err),
mpredVar_err = mean(predVar_err),
mpredRP_rerr = mean(predRP_rerr),
mpredVar_rerr = mean(predVar_rerr),
cv= mean(cv),
pred_cv = mean(pred_cv),
mpredcv_err = mean(predcv_err),
mpredcv_rerr = mean(predcv_rerr))
m_predY_acc = mpredY_sub%>%
dplyr::group_by(Exp) %>%
dplyr::summarize(mpred_rerr = mean(predRP_rerr)*100,
mpredVar_rerr = mean(predVar_rerr)*100,
mpredcv_rerr = mean(predcv_rerr)*100)
m_predY_acc
WAIC and LOO-CV
#extract waic and loo-cv from parameter list
m_WAIC <- dplyr::group_by(Bayparlist, Exp) %>%
dplyr::summarize(n =n(),
m_looic = mean(looic),
m_waic = mean(waic),
se_waic = sd(waic)/sqrt(n-1),
se_looic = sd(looic)/sqrt(n-1),
m_p_loo = mean(p_loo),
m_elpd_loo = mean(elpd_loo),
m_se_looic = mean(se_looic),
m_se_p_loo = mean(se_p_loo),
m_p_waic = mean(p_waic),
m_se_waic = mean(se_waic))
m_WAIC
#load test results
AllDat_newY <- read.csv(paste0(modelPath, "/rlt/AllDat_newY_",modelversion,".csv"))
AllDat_newY$WMSize <- as.factor(AllDat_newY$WMSize)
levels(AllDat_newY$WMSize) = c("low", "medium", "high")
AllDat_newY[which(AllDat_newY$Exp == "Exp4"), "Exp"] = "Exp4a"
AllDat_newY[which(AllDat_newY$Exp == "Exp5"), "Exp"] = "Exp4b"
Plot Results
RP biase
RP_bias <- ggplot(data = m_predY, aes(x = curDur, y = m_m_repDur - curDur,color=WMSize, shape = as.factor('Observation'))) +
geom_point(size=2, alpha = 0.5)+
geom_line(data= m_newY, aes(x=curDur, y=m_mu_r-curDur, color=WMSize)) +
geom_hline(yintercept = 0, linetype='dashed')+
facet_grid(cols = vars(Exp)) +
labs(x=" ", y="Reproduction bias (s)", shape=" ", color = "Memory Load")+theme_new+
colorSet3+guides(shape="none")
ggsave(paste0(getwd(), "/", modelPath, "/figures/RP_bias.png"), RP_bias, width = 6, height = 6)
RP_bias
RP_bias_obs <- ggplot(data = AllExpData%>% filter(Exp != 'Exp5') %>%
dplyr::group_by(Exp, curDur, WMSize, NSub) %>%
dplyr::summarize(n = n(),
m_repDur = mean(repDur),
se_repDur = sd(repDur)/sqrt(n-1)) %>%
dplyr::group_by(Exp, curDur, WMSize) %>%
dplyr::summarize(m_m_repDur = mean(m_repDur),
m_se_repDur = mean(se_repDur)), aes(x = curDur, y = m_m_repDur-curDur, color=as.factor(WMSize))) +
geom_point()+
geom_line()+
#geom_errorbar(width=.2, aes(ymin = m_m_repDur-curDur - m_se_repDur, ymax = m_m_repDur -curDur + m_se_repDur)) +
geom_hline(yintercept = 0, linetype='dashed')+
facet_grid(cols = vars(Exp)) +
labs(x="Sample intervals (s)", y="Reproduction bias(s)", shape=" ", color = "Memory Load")+
theme_new+colorSet3+guides(shape="none")+
scale_x_continuous(breaks=seq(0, 1.6, 0.4))+ theme(legend.position="top")
## `summarise()` has grouped output by 'Exp', 'curDur', 'WMSize'. You can override using the `.groups` argument.
## `summarise()` has grouped output by 'Exp', 'curDur'. You can override using the `.groups` argument.
ggsave(paste0(getwd(), "/", modelPath, "/figures/RP_bias_obs.png"), RP_bias_obs, width = 6, height = 4)
RP_bias_obs
RP_bias_obs <- ggplot(data = AllExpData%>% filter(Exp != 'Exp5') %>%
dplyr::group_by(Exp, curDur, WMSize, NSub) %>%
dplyr::summarize(n = n(),
m_repDur = mean(repDur),
se_repDur = sd(repDur)/sqrt(n-1)) %>%
dplyr::group_by(Exp, curDur, WMSize) %>%
dplyr::summarize(m_m_repDur = mean(m_repDur),
m_se_repDur = mean(se_repDur)), aes(x = curDur, y = m_m_repDur-curDur, color=as.factor(WMSize))) +
geom_point()+
geom_line()+
geom_errorbar(width=.2, aes(ymin = m_m_repDur-curDur - m_se_repDur, ymax = m_m_repDur -curDur + m_se_repDur)) +
geom_hline(yintercept = 0, linetype='dashed')+
facet_grid(cols = vars(Exp)) +
labs(x="Sample intervals (s)", y="Reproduction bias(s)", shape=" ", color = "Memory Load")+
theme_new+colorSet3+guides(shape="none")+
scale_x_continuous(breaks=seq(0, 1.6, 0.4))+ theme(legend.position="top")
## `summarise()` has grouped output by 'Exp', 'curDur', 'WMSize'. You can override using the `.groups` argument.
## `summarise()` has grouped output by 'Exp', 'curDur'. You can override using the `.groups` argument.
ggsave(paste0(getwd(), "/", modelPath, "/figures/RP_bias_obs.png"), RP_bias_obs, width = 6, height = 4)
RP_bias_obs
RP <- ggplot(data = m_predY, aes(x = curDur, y = m_m_repDur, color=WMSize, shape = as.factor('Observation'))) +
geom_point(size=2, alpha = 0.5)+
geom_line(data= m_newY, aes(x=curDur, y=m_mu_r, color=WMSize)) +
geom_abline(slope=1, intercept=0)+
facet_grid(cols = vars(Exp)) +
labs(x="Sample intervals (s)", y="Reproduction (s)", shape=" ", color = "Memory Load")+
theme_new+colorSet3+guides(shape="none")+theme_new +theme(legend.position="top")
ggsave(paste0(getwd(), "/", modelPath, "/figures/RP.png"), RP, width = 6, height = 6)
RP
RP Slope
# calculate the slope of the cv curve
RPslope_model <- function(df) {
lm(m_repDur ~ log(curDur), data = df)
}
RPslopes <- mpredY_sub %>%
dplyr::group_by(NSub, Exp, WMSize) %>% nest() %>%
mutate(model = map(data, RPslope_model)) %>%
mutate(slope = map(model, broom::tidy)) %>% # get estimates
unnest(slope, .drop = TRUE) %>% # remove raw data
select(-std.error,-statistic, -p.value) %>% # remove unnessary columns
spread(term, estimate) %>% # spread stimates
dplyr::rename(Intercept = `(Intercept)`, slope = `log(curDur)`)
RPslopes$data <- NULL
RPslopes$model <- NULL
mRPslopes <- RPslopes%>% dplyr::group_by(WMSize, Exp) %>%
dplyr::summarize(m_Intercept = mean(Intercept),
m_slope = mean(slope),
n = n(),
se_slope = sd(slope)/sqrt(n-1),
se_Intercept = sd(Intercept)/sqrt(n-1))
## `summarise()` has grouped output by 'WMSize'. You can override using the
## `.groups` argument.
plt_CVslope <- ggplot(mRPslopes, aes(Exp, m_slope, ymin = m_slope - se_slope, ymax = m_slope + se_slope, group =interaction(Exp, WMSize), color = factor(WMSize)), shape = factor(WMSize))+
geom_line(stat = "identity",position = position_dodge(width = 0.2))+
geom_point(stat = "identity",position = position_dodge(width = 0.2))+
geom_errorbar(width=.2, position = position_dodge(width = 0.2)) +
colorSet5+
labs(x = "", y = TeX("Slope of RP"), color = 'Memory Load') +
theme_new +theme(legend.position="top")
plt_CVslope
## geom_path: Each group consists of only one observation. Do you need to adjust
## the group aesthetic?
Anova analysis on slopes of RP
RPslopes$WMSize = as.factor(RPslopes$WMSize)
ezANOVA(data = RPslopes, dv= slope, wid=NSub, within=.(WMSize), between = .(Exp))
## Warning: Converting "NSub" to factor for ANOVA.
## Warning: Converting "Exp" to factor for ANOVA.
## Warning: The column supplied as the wid variable contains non-unique values
## across levels of the supplied between-Ss variables. Automatically fixing this by
## generating unique wid labels.
## $ANOVA
## Effect DFn DFd F p p<.05 ges
## 2 Exp 4 75 2.465821 5.210071e-02 0.11127986
## 3 WMSize 2 150 24.943347 4.446854e-10 * 0.01567449
## 4 Exp:WMSize 8 150 6.439500 3.455596e-07 * 0.01617814
##
## $`Mauchly's Test for Sphericity`
## Effect W p p<.05
## 3 WMSize 0.943444 0.1160104
## 4 Exp:WMSize 0.943444 0.1160104
##
## $`Sphericity Corrections`
## Effect GGe p[GG] p[GG]<.05 HFe p[HF] p[HF]<.05
## 3 WMSize 0.9464714 1.164847e-09 * 0.9702459 7.594066e-10 *
## 4 Exp:WMSize 0.9464714 6.643447e-07 * 0.9702459 4.968822e-07 *
ezANOVA(data = RPslopes %>% filter(Exp =='Exp1'), dv= slope, wid=NSub, within = .(WMSize))
## Warning: Converting "NSub" to factor for ANOVA.
## $ANOVA
## Effect DFn DFd F p p<.05 ges
## 2 WMSize 2 30 1.628257 0.2131416 0.002398051
##
## $`Mauchly's Test for Sphericity`
## Effect W p p<.05
## 2 WMSize 0.8320858 0.2761702
##
## $`Sphericity Corrections`
## Effect GGe p[GG] p[GG]<.05 HFe p[HF] p[HF]<.05
## 2 WMSize 0.8562273 0.2172718 0.9557549 0.214484
ezANOVA(data = RPslopes %>% filter(Exp =='Exp2'), dv= slope, wid=NSub, within = .(WMSize))
## Warning: Converting "NSub" to factor for ANOVA.
## $ANOVA
## Effect DFn DFd F p p<.05 ges
## 2 WMSize 2 30 21.20007 1.822755e-06 * 0.07646424
##
## $`Mauchly's Test for Sphericity`
## Effect W p p<.05
## 2 WMSize 0.8461525 0.3105563
##
## $`Sphericity Corrections`
## Effect GGe p[GG] p[GG]<.05 HFe p[HF] p[HF]<.05
## 2 WMSize 0.8666656 7.302883e-06 * 0.9698478 2.493639e-06 *
ezANOVA(data = RPslopes %>% filter(Exp =='Exp3'), dv= slope, wid=NSub, within = .(WMSize))
## Warning: Converting "NSub" to factor for ANOVA.
## $ANOVA
## Effect DFn DFd F p p<.05 ges
## 2 WMSize 2 30 0.4137346 0.6648914 0.00152841
##
## $`Mauchly's Test for Sphericity`
## Effect W p p<.05
## 2 WMSize 0.9800567 0.8684769
##
## $`Sphericity Corrections`
## Effect GGe p[GG] p[GG]<.05 HFe p[HF] p[HF]<.05
## 2 WMSize 0.9804466 0.6610089 1.126392 0.6648914
ezANOVA(data = RPslopes %>% filter(Exp =='Exp4a'), dv= slope, wid=NSub, within = .(WMSize))
## Warning: Converting "NSub" to factor for ANOVA.
## $ANOVA
## Effect DFn DFd F p p<.05 ges
## 2 WMSize 2 30 7.246942 0.002705914 * 0.01861911
##
## $`Mauchly's Test for Sphericity`
## Effect W p p<.05
## 2 WMSize 0.9730381 0.8258648
##
## $`Sphericity Corrections`
## Effect GGe p[GG] p[GG]<.05 HFe p[HF] p[HF]<.05
## 2 WMSize 0.973746 0.002974256 * 1.117022 0.002705914 *
ezANOVA(data = RPslopes %>% filter(Exp =='Exp4b'), dv= slope, wid=NSub, within = .(WMSize))
## Warning: Converting "NSub" to factor for ANOVA.
## $ANOVA
## Effect DFn DFd F p p<.05 ges
## 2 WMSize 2 30 8.551317 0.001151174 * 0.03802106
##
## $`Mauchly's Test for Sphericity`
## Effect W p p<.05
## 2 WMSize 0.8930404 0.4529994
##
## $`Sphericity Corrections`
## Effect GGe p[GG] p[GG]<.05 HFe p[HF] p[HF]<.05
## 2 WMSize 0.9033753 0.001756187 * 1.019764 0.001151174 *
Anova analysis on Intercept of RP
ezANOVA(data = RPslopes, dv= Intercept, wid=NSub, within=.(WMSize), between = .(Exp))
## Warning: Converting "NSub" to factor for ANOVA.
## Warning: Converting "Exp" to factor for ANOVA.
## Warning: The column supplied as the wid variable contains non-unique values
## across levels of the supplied between-Ss variables. Automatically fixing this by
## generating unique wid labels.
## $ANOVA
## Effect DFn DFd F p p<.05 ges
## 2 Exp 4 75 0.6157298 6.526476e-01 0.028876673
## 3 WMSize 2 150 4.6236001 1.125794e-02 * 0.005792574
## 4 Exp:WMSize 8 150 6.6943967 1.762289e-07 * 0.032641724
##
## $`Mauchly's Test for Sphericity`
## Effect W p p<.05
## 3 WMSize 0.8020579 0.0002855036 *
## 4 Exp:WMSize 0.8020579 0.0002855036 *
##
## $`Sphericity Corrections`
## Effect GGe p[GG] p[GG]<.05 HFe p[HF] p[HF]<.05
## 3 WMSize 0.8347649 1.634417e-02 * 0.8515168 1.573692e-02 *
## 4 Exp:WMSize 0.8347649 1.479632e-06 * 0.8515168 1.191853e-06 *
CV
curDurItem <- unique(m_predY$curDur)
RP_CV <- ggplot(data= m_predY, aes(x=curDur, y= cv, color=WMSize, shape = as.factor('Observation'))) +
geom_point(size=2, alpha = 0.5)+
geom_line(data = m_newY, aes(x=curDur, y= m_sig_r/m_mu_r, color=WMSize)) +
facet_grid(~Exp) +
labs(x="Duration (s)", y="CV", shape=" ", color = "Memory Load")+ theme_new+
colorSet3+guides(shape="none")+theme(strip.text.x = element_blank())
RP_CV
ggsave(paste0(getwd(), "/", modelPath, "/figures/RP_CV.png"), RP_CV, width = 6, height = 6)
CV slope
# calculate the slope of the cv curve
cvSlope_model <- function(df) {
lm(log(cv_obs) ~ log(curDur), data = df)
}
mpredY <- dplyr::group_by(AllDat_predY, curDur, WMSize, Exp, NSub) %>%
dplyr::summarize(m_repDur = mean(repDur),
sd_repDur = sd(repDur),
n = n(),
sd_repDur = sd(repDur),
m_mu_r = mean(mu_r),
m_sig_r = mean(sig_r),
wp = mean(wp),
log_lik =mean(log_lik))
## `summarise()` has grouped output by 'curDur', 'WMSize', 'Exp'. You can override
## using the `.groups` argument.
mpredY$cv_obs <- mpredY$sd_repDur/mpredY$m_repDur
CVslopes <- mpredY %>%
dplyr::group_by(NSub, Exp, WMSize) %>% nest() %>% # nested data
mutate(model = map(data, cvSlope_model)) %>% # linear regression
mutate(slope = map(model, broom::tidy)) %>% # get estimates
unnest(slope, .drop = TRUE) %>% # remove raw data
select(-std.error,-statistic, -p.value) %>% # remove unnessary columns
spread(term, estimate) %>% # spread stimates
dplyr::rename(obs_cv_Intercept = `(Intercept)`, obs_cv_slope = `log(curDur)`) # rename columns
CVslopes$data = NULL
CVslopes$model = NULL
#change the table struction of slopes for spss
CVslopes_list <-split(CVslopes, CVslopes$WMSize)
CVslopes_spss = NULL
for (i in 1: length(CVslopes_list)){
temp = CVslopes_list[[i]]
WMSize = unique(temp$WMSize)
temp$WMSize = NULL
colnames(temp) = c('Exp', 'NSub', paste0('obs_cv_Intercept_',WMSize), paste0('obs_cv_slope_',WMSize))
if(i == 1)
CVslopes_spss = temp
else
CVslopes_spss = left_join(CVslopes_spss, temp, by=c("Exp", "NSub"))
}
mCVslopes <- CVslopes%>% dplyr::group_by(WMSize, Exp) %>%
dplyr::summarize(n =n(),
m_obs_cv_Intercept = mean(obs_cv_Intercept),
m_obs_cv_slope = mean(obs_cv_slope),
se_obs_CV_slope = sd(obs_cv_slope)/sqrt(n-1),
se_obs_CV_Intercept = sd(obs_cv_Intercept)/sqrt(n-1))
## `summarise()` has grouped output by 'WMSize'. You can override using the
## `.groups` argument.
plt_CVslope <- ggplot(mCVslopes, aes(Exp, m_obs_cv_slope, ymin = m_obs_cv_slope - se_obs_CV_slope, ymax = m_obs_cv_slope + se_obs_CV_slope, group =interaction(Exp, WMSize), color = factor(WMSize)), shape = factor(WMSize))+
geom_line(stat = "identity",position = position_dodge(width = 0.2))+
geom_point(stat = "identity",position = position_dodge(width = 0.2))+
geom_errorbar(width=.2, position = position_dodge(width = 0.2)) +
colorSet5+
labs(x = "", y = TeX("Slope of CV"), color = 'Memory Load') +
theme_new +theme(legend.position="top")
plt_CVslope
## geom_path: Each group consists of only one observation. Do you need to adjust
## the group aesthetic?
ggsave(paste0(getwd(), "/", modelPath, "/figures/plt_CVslope.png"), plt_CVslope, width = 4, height = 4)
## geom_path: Each group consists of only one observation. Do you need to adjust
## the group aesthetic?
plt_CVIntecept<- ggplot(mCVslopes, aes(Exp, m_obs_cv_Intercept, ymin = m_obs_cv_Intercept - se_obs_CV_Intercept, ymax = m_obs_cv_Intercept + se_obs_CV_Intercept, group =interaction(Exp, WMSize), color = factor(WMSize)), shape = factor(WMSize))+
geom_line(stat = "identity",position = position_dodge(width = 0.2))+
geom_point(stat = "identity",position = position_dodge(width = 0.2))+
geom_errorbar(width=.2, position = position_dodge(width = 0.2)) +
colorSet5+
labs(x = "", y = TeX("Intercept of CV"), color = 'Memory Load') +
theme_new +theme(legend.position="top")
plt_CVIntecept
## geom_path: Each group consists of only one observation. Do you need to adjust
## the group aesthetic?
ggsave(paste0(getwd(), "/", modelPath, "/figures/plt_Intercept.png"), plt_CVIntecept, width = 4, height = 4)
## geom_path: Each group consists of only one observation. Do you need to adjust
## the group aesthetic?
## Figures in the MS
fig_CV_slope_Intecept<-ggarrange(plt_CVslope, plt_CVIntecept, common.legend = TRUE, ncol=2, nrow=1, labels = c("a", "b"))
## geom_path: Each group consists of only one observation. Do you need to adjust
## the group aesthetic?
## geom_path: Each group consists of only one observation. Do you need to adjust
## the group aesthetic?
## geom_path: Each group consists of only one observation. Do you need to adjust
## the group aesthetic?
ggsave(paste0(getwd(), "/", modelPath, "/figures/fig_CV_slope_Intecept.png"), fig_CV_slope_Intecept, width = 8, height = 4)
fig_CV_slope_Intecept
## Figures in the MS
fig3<-ggarrange(RP_bias, RP_CV, common.legend = TRUE, ncol=1, nrow=2, labels = c("a", "b"))
ggsave(paste0(getwd(), "/", modelPath, "/figures/fig3.png"), fig3, width = 6, height = 5)
fig3
Indifference Pointtemp_data = sumexpSub_diff %>%
Indifference Point(linear regression)
#Observed Indifference Point
obs_model <- function(df) {
lm(repDur ~ curDur, data = df)
}
obs_InP <- AllExpData %>% #filter(curDur %in% c(1.1,1.4)) %>%
dplyr::group_by(NSub, Exp, WMSize) %>% nest() %>%
mutate(model = map(data, obs_model)) %>% # linear regression
mutate(slope = map(model, broom::tidy)) %>% # get estimates
unnest(slope, .drop = TRUE) %>% # remove raw data
select(-std.error,-statistic, -p.value) %>% # remove unnessary columns
spread(term, estimate) %>% # spread stimates
dplyr::rename(Intercept = `(Intercept)`, slope = curDur) # rename columns
obs_InP$model = NULL
obs_InP$data = NULL
obs_InP$inP = obs_InP$Intercept /(1-obs_InP$slope)
write.csv(obs_InP, file = paste0(modelPath, "/rlt/obs_InP.csv"))
mobs_InP = obs_InP %>% group_by(Exp, WMSize)%>%
dplyr::summarise(n=n(),
m_Intercept = mean(Intercept),
se_Intercept= sd(Intercept)/sqrt(n-1),
m_inP = mean(inP),
se_inP = sd(inP)/sqrt(n-1),
m_slope = mean(slope),
se_slope = sd(slope)/sqrt(n-1))
## `summarise()` has grouped output by 'Exp'. You can override using the `.groups`
## argument.
plt_InP_linear<- ggplot(data = mobs_InP, aes(x= Exp, y=m_inP, color = WMSize))+
geom_line(stat = "identity",position = position_dodge(width = 0.2))+
geom_point(stat = "identity",position = position_dodge(width = 0.2))+
geom_errorbar(width=.3, aes(ymin = m_inP - se_inP, ymax = m_inP + se_inP), position = position_dodge(width = 0.2)) +theme_new+
labs(colour = "Memory Load")+colorSet3+
xlab(' ')+ylab("indifference point (s)")+guides(shape="none")+
theme(legend.position = "top")
ggsave(paste0(getwd(), "/", modelPath, "/figures/plt_InP_linear.png"), plt_InP_linear, width = 5, height = 5)
## geom_path: Each group consists of only one observation. Do you need to adjust
## the group aesthetic?
plt_InP_linear
## geom_path: Each group consists of only one observation. Do you need to adjust
## the group aesthetic?
RP_Exp2 <- ggplot(data = m_predY%>% filter(Exp == 'Exp2'), aes(x = curDur, y = m_m_repDur, color=WMSize, shape = as.factor('Observation'))) +
geom_point(size=2, alpha = 0.5)+
geom_segment(data =mobs_InP%>%filter(Exp == 'Exp2'), aes(x = 0.3, y = 0.3*m_slope+m_Intercept, xend = 2.2, yend = 2.2*m_slope+m_Intercept), arrow = NULL)+ ##dotted lines
geom_point(data =mobs_InP%>% filter(Exp == 'Exp2'), aes(x = m_inP, y = 0.1, color=WMSize), shape=21) + ## Intersection
geom_segment(data =mobs_InP%>% filter(Exp == 'Exp2'), aes(x=m_inP, y = 0.1, xend = m_inP, yend = m_inP), arrow = NULL)+ ##vertical lines for Intersection
#geom_line(data= m_newY%>% filter(Exp == 'Exp2'), aes(x=curDur, y=m_mu_r, color=WMSize)) +
geom_abline(slope=1, intercept=0, linetype = 2)+
labs(x="Sample intervals (s)", y="Reproduction (s)", shape=" ", color = "Memory Load")+
theme_new+colorSet3+guides(shape="none")+ ylim(0,2)
ggsave(paste0(getwd(), "/", modelPath, "/figures/RP_Exp2.png"), RP_Exp2, width = 4, height = 4)
RP_Exp2
Indifference Point (bootstraps)
#### Estimated Indifference Point
mRep_model <- function(df) {
lm(mu_r ~ curDur, data = df)
}
#Observed Indifference Point
obs_model <- function(df) {
lm(repDur ~ curDur, data = df)
}
# Custom function to find predicted indifference point
getPredInP_boot <- function(df, idx){
vars <- c('NSub', 'Exp', 'WMSize')
gp_vars = syms(vars)
slopes <- df[idx, ] %>%
dplyr::group_by(!!!gp_vars) %>% nest() %>% # nested data
mutate(model = map(data, mRep_model)) %>% # linear regression
mutate(slope = map(model, broom::tidy)) %>% # get estimates out
unnest(slope, .drop = TRUE) %>% # remove raw data
select(-std.error,-statistic, -p.value) %>% # remove unnessary clumns
spread(term, estimate) %>% # spread stimates
dplyr::rename(minRP = `(Intercept)`, slope = curDur) # rename columns
slopes$inP = slopes$minRP /(1-slopes$slope)
return(slopes$inP)
}
# Custom function to find observed indifference point
getRPInP_boot <- function(df, idx){
vars <- c('NSub', 'Exp', 'WMSize')
gp_vars = syms(vars)
slopes <- df[idx, ] %>%
dplyr::group_by(!!!gp_vars) %>% nest() %>% # nested data
mutate(model = map(data, obs_model)) %>% # linear regression
mutate(slope = map(model, broom::tidy)) %>% # get estimates out
unnest(slope, .drop = TRUE) %>% # remove raw data
select(-std.error,-statistic, -p.value) %>% # remove unnessary clumns
spread(term, estimate) %>% # spread stimates
dplyr::rename(minRP = `(Intercept)`, slope = curDur) # rename columns
slopes$inP = slopes$minRP /(1-slopes$slope)
return(slopes$inP)
}
#calculate the bootstrapped 95% confidence intervals
generateCI = FALSE # tag for generation CI
if(generateCI){
cilist <- NULL
for(expname in unique(AllDat_predY$Exp)){
for(nsub in unique(AllDat_predY$NSub)){
for(WMsize in unique(AllDat_predY$WMSize)){
dat = AllDat_predY %>% filter(Exp == expname, NSub == nsub, WMSize ==WMsize)
set.seed(100)
num = 1000
bs_predInP <- boot(dat, getPredInP_boot, R = num)
bs_RPInP <- boot(dat, getRPInP_boot, R = num)
ci = data.frame(
Exp = expname,
NSub = nsub,
WMSize = WMsize,
sd_predInP_boot = sd(bs_predInP$t),
m_predInP_boot = median(bs_predInP$t),
sd_RPInP_boot = sd(bs_RPInP$t),
mRP_InP_boot = median(bs_RPInP$t)
)
cilist = data.frame(rbind(cilist, ci))
}
}
}
write.csv(cilist, file = paste0("ci_list_median_1000.csv"))
}
# load the generated indifference point values and mark the outlier
cilist = read.csv(paste0("ci_list_median_1000.csv"))
cilist$Exp = as.factor(cilist$Exp)
cilist$WMSize= as.factor(cilist$WMSize)
cilist$inPOutlier = FALSE
cilist[which(cilist$mRP_InP_boot > 1.7 |cilist$mRP_InP_boot < 0.5 | cilist$m_predInP_boot < 0.5|cilist$m_predInP_boot > 1.7),"inPOutlier"] = TRUE
#check if the outlier is the same as the outliers in variable slope_pr
cilist%>% filter(inPOutlier == TRUE)
Inifference Point (curve)
AllDat_newY$predErr = AllDat_newY$mu_r - AllDat_newY$curDur
temp_newY <- AllDat_newY %>% filter(curDur > 0.8, curDur < 1.1) %>% select(Exp, WMSize, NSub, predErr, curDur)
InP_curve<- temp_newY%>% dplyr::group_by(Exp, WMSize, NSub)%>%
dplyr::summarise(minErr = min(abs(predErr)), idx = which.min(abs(predErr)))
## `summarise()` has grouped output by 'Exp', 'WMSize'. You can override using the
## `.groups` argument.
InP_curve$InP_curve = temp_newY[InP_curve$idx,]$curDur
InP_curve$y = temp_newY[InP_curve$idx,]$predErr + temp_newY[InP_curve$idx,]$curDur
InP_curve
Inifference Point (Equation)
#sig_sm2 = sig_s2 + ls*(size[i]);
#wp_pr[i] = sig_sm2 /(sig_sm2 + sigma_pr2);
Bayparlist$wp_1 = (Bayparlist$sig_s2 + Bayparlist$ls *1 )/ (Bayparlist$sig_s2 + Bayparlist$ls *1 + Bayparlist$sig_pr2_log)
Bayparlist$wp_3 = (Bayparlist$sig_s2 + Bayparlist$ls *2) / (Bayparlist$sig_s2 + Bayparlist$ls *2 + Bayparlist$sig_pr2_log)
Bayparlist$wp_5 = (Bayparlist$sig_s2 + Bayparlist$ls *3) / (Bayparlist$sig_s2 + Bayparlist$ls *3 + Bayparlist$sig_pr2_log)
#log(IP)=[kr*M + (1-w_p)ks*M+ w_p*mu_p]/w_p
size = 1
Bayparlist$InP_1 = exp((Bayparlist$kr*size + (1-Bayparlist$wp_1)*Bayparlist$ks*size + Bayparlist$wp_1* Bayparlist$mu_pr_log)/Bayparlist$wp_1)
size = 2
Bayparlist$InP_3 = exp((Bayparlist$kr*size + (1-Bayparlist$wp_3)*Bayparlist$ks*size + Bayparlist$wp_3* Bayparlist$mu_pr_log)/Bayparlist$wp_3)
size = 3
Bayparlist$InP_5 = exp((Bayparlist$kr*size + (1-Bayparlist$wp_5)*Bayparlist$ks*size + Bayparlist$wp_5* Bayparlist$mu_pr_log)/Bayparlist$wp_5)
Bayparlist$sig2_post_1 = (Bayparlist$sig_pr2_log *(Bayparlist$sig_s2 +Bayparlist$ls*1))/(Bayparlist$sig_pr2_log + Bayparlist$sig_s2 +Bayparlist$ls*1) *0.5
Bayparlist$sig2_post_3 =(Bayparlist$sig_pr2_log *(Bayparlist$sig_s2 +Bayparlist$ls*3))/(Bayparlist$sig_pr2_log + Bayparlist$sig_s2 +Bayparlist$ls*3) *0.5
Bayparlist$sig2_post_5 = (Bayparlist$sig_pr2_log *(Bayparlist$sig_s2 +Bayparlist$ls*5))/(Bayparlist$sig_pr2_log + Bayparlist$sig_s2 +Bayparlist$ls*5) *0.5
Bayparlist%>%filter(Exp == 'Exp2')%>%select("NSub", "Exp", "sig2_post_1", "sig2_post_3","sig2_post_5")
Bayparlist$InP_old_1 = exp((Bayparlist$kr*size - (1-Bayparlist$wp_1)*Bayparlist$ks*size +(Bayparlist$sig_pr2_log *(Bayparlist$sig_s2 +Bayparlist$ls*1))/(Bayparlist$sig_pr2_log + Bayparlist$sig_s2 +Bayparlist$ls*1) *0.5 + Bayparlist$wp_1* Bayparlist$mu_pr_log)/Bayparlist$wp_1)
Bayparlist$InP_old_3 = exp((Bayparlist$kr*size - (1-Bayparlist$wp_3)*Bayparlist$ks*size +(Bayparlist$sig_pr2_log *(Bayparlist$sig_s2 +Bayparlist$ls*3))/(Bayparlist$sig_pr2_log + Bayparlist$sig_s2 +Bayparlist$ls*3) *0.5 + Bayparlist$wp_3* Bayparlist$mu_pr_log)/Bayparlist$wp_3)
Bayparlist$InP_old_5 = exp((Bayparlist$kr*size - (1-Bayparlist$wp_5)*Bayparlist$ks*size +(Bayparlist$sig_pr2_log *(Bayparlist$sig_s2 +Bayparlist$ls*5))/(Bayparlist$sig_pr2_log + Bayparlist$sig_s2 +Bayparlist$ls*5) *0.5 + Bayparlist$wp_5* Bayparlist$mu_pr_log)/Bayparlist$wp_5)
Bayparlist%>%select(NSub, Exp, wp_1, wp_3, wp_5, InP_1, InP_3, InP_5, InP_old_1, InP_old_3, InP_old_5)%>%
unite(newcol1, wp_1, InP_1, InP_old_1)%>%
unite(newcol3, wp_3, InP_3, InP_old_3)%>%
unite(newcol5, wp_5, InP_5, InP_old_5)%>%
melt(id.vars = c("NSub", "Exp")) %>%
dplyr::rename(
WMSize = variable,
newcol = value
)%>%
separate(newcol, c("wp", "InP_Eq", "InP_Eq_old"), sep = "_") -> Baypar_WMSize
Baypar_WMSize$WMSize <- substr(Baypar_WMSize$WMSize, 7, 7)
Baypar_WMSize$WMSize <- as.factor(Baypar_WMSize$WMSize)
Baypar_WMSize$InP_Eq = as.numeric(Baypar_WMSize$InP_Eq)
Baypar_WMSize$InP_Eq_old = as.numeric(Baypar_WMSize$InP_Eq_old)
levels(Baypar_WMSize$WMSize) = c("low", "medium", "high")
##save csv for SPSS
left_join(Baypar_WMSize, InP_curve%>%select("Exp", "WMSize","NSub", "InP_curve"), by = c("Exp","WMSize","NSub")) -> Baypar_WMSize
left_join(Baypar_WMSize, cilist %>% select("Exp", "NSub", "WMSize", "m_predInP_boot", "mRP_InP_boot"), by = c("Exp","WMSize","NSub"))-> Baypar_WMSize
check indifference Points
Baypar_WMSize_melt <- melt(Baypar_WMSize%>%select("Exp","WMSize","NSub", "InP_Eq", "InP_Eq_old", "InP_curve", "m_predInP_boot", "mRP_InP_boot"), id.vars = c("Exp","WMSize","NSub"),
variable.name = "Type",
value.name = "InP" )
Baypar_WMSize_melt %>%filter(Type== "InP_curve") %>% dplyr::group_by(Exp)%>%
dplyr::summarise(n = n(),
m_InP = mean(InP), se_InP= sd(InP)/sqrt(n-1))
Baypar_WMSize_melt$Exp = as.factor(Baypar_WMSize_melt$Exp)
Baypar_WMSize_melt%>% dplyr::group_by(Exp, WMSize, Type)%>%
dplyr::summarise(n = n(),
m_InP = mean(InP), se_InP= sd(InP)/sqrt(n-1))-> mBaypar_InP
## `summarise()` has grouped output by 'Exp', 'WMSize'. You can override using the
## `.groups` argument.
#plot InP_curve(the intersections of the Prediction curve with the diagonal)
plt_InP<- ggplot(data = mBaypar_InP%>%filter(Type== "InP_curve"), aes(x= Exp, y=m_InP, color = WMSize, shape = Type))+
geom_line(stat = "identity",position = position_dodge(width = 0.2))+
geom_point(stat = "identity",position = position_dodge(width = 0.2))+
geom_errorbar(width=.2, aes(ymin = m_InP - se_InP, ymax = m_InP + se_InP), position = position_dodge(width = 0.2)) +theme_new+
labs(colour = "Memory Load")+colorSet3+
xlab(' ')+ylab("indifference point (s)")+guides(shape="none")+
theme(legend.position = "top")
ggsave(paste0(getwd(), "/", modelPath, "/figures/plt_InP.png"), plt_InP, width = 3, height = 3)
## geom_path: Each group consists of only one observation. Do you need to adjust
## the group aesthetic?
plt_InP
## geom_path: Each group consists of only one observation. Do you need to adjust
## the group aesthetic?
temp_inp = mBaypar_InP%>%filter(Type== "InP_curve", Exp == 'Exp2')
RP_bias_Exp2 <- ggplot(data = m_predY%>%filter(Exp == 'Exp2'), aes(x = curDur, y = m_m_repDur - curDur,color=WMSize, shape = as.factor('Observation'))) +
geom_point(size=2, alpha = 0.5)+
geom_line(data= m_newY%>%filter(Exp == 'Exp2'), aes(x=curDur, y=m_mu_r-curDur, color=WMSize)) +
geom_hline(yintercept = 0, linetype='dashed')+
geom_point(data =temp_inp, aes(x = m_InP+0.03, y = -0.5, color=WMSize), shape=21) + ## Intersection
geom_segment(data =temp_inp, aes(x=m_InP+0.03, y = -0.5, xend = m_InP+0.03, yend = 0), arrow = NULL)+ ##vertical lines for Intersection
labs(x=" ", y="Reproduction bias (s)", shape=" ", color = "Memory Load")+theme_new+
colorSet3+guides(shape="none")
ggsave(paste0(getwd(), "/", modelPath, "/figures/RP_bias_Exp2.png"), RP_bias_Exp2, width = 5, height = 3)
RP_bias_Exp2
anova on observed InP
ezANOVA(data = Baypar_WMSize, dv= InP_curve, wid=NSub, within = .(WMSize), between = .(Exp))
## Warning: Converting "NSub" to factor for ANOVA.
## Warning: Converting "Exp" to factor for ANOVA.
## Warning: The column supplied as the wid variable contains non-unique values
## across levels of the supplied between-Ss variables. Automatically fixing this by
## generating unique wid labels.
## $ANOVA
## Effect DFn DFd F p p<.05 ges
## 2 Exp 4 75 1.122790 3.522214e-01 0.05010528
## 3 WMSize 2 150 4.862201 8.994962e-03 * 0.00766407
## 4 Exp:WMSize 8 150 5.529029 3.939989e-06 * 0.03393766
##
## $`Mauchly's Test for Sphericity`
## Effect W p p<.05
## 3 WMSize 0.1284775 1.063161e-33 *
## 4 Exp:WMSize 0.1284775 1.063161e-33 *
##
## $`Sphericity Corrections`
## Effect GGe p[GG] p[GG]<.05 HFe p[HF] p[HF]<.05
## 3 WMSize 0.5343243 0.0280800589 * 0.535749 0.0279830198 *
## 4 Exp:WMSize 0.5343243 0.0004160714 * 0.535749 0.0004100909 *
pairwise.t.test(Baypar_WMSize$InP_curve, Baypar_WMSize$Exp)
##
## Pairwise comparisons using t tests with pooled SD
##
## data: Baypar_WMSize$InP_curve and Baypar_WMSize$Exp
##
## Exp1 Exp2 Exp3 Exp4a
## Exp2 1.000 - - -
## Exp3 1.000 1.000 - -
## Exp4a 1.000 1.000 1.000 -
## Exp4b 0.125 0.030 0.049 0.690
##
## P value adjustment method: holm
ezANOVA(data = Baypar_WMSize %>% filter(Exp =='Exp2'), dv= InP_curve, wid=NSub, within = .(WMSize))
## Warning: Converting "NSub" to factor for ANOVA.
## $ANOVA
## Effect DFn DFd F p p<.05 ges
## 2 WMSize 2 30 6.278925 0.005273799 * 0.04191182
##
## $`Mauchly's Test for Sphericity`
## Effect W p p<.05
## 2 WMSize 0.1557357 2.221859e-06 *
##
## $`Sphericity Corrections`
## Effect GGe p[GG] p[GG]<.05 HFe p[HF] p[HF]<.05
## 2 WMSize 0.5422216 0.02127963 * 0.5515802 0.02067639 *
ezANOVA(data = Baypar_WMSize %>% filter(Exp =='Exp3'), dv= InP_curve, wid=NSub, within = .(WMSize))
## Warning: Converting "NSub" to factor for ANOVA.
## $ANOVA
## Effect DFn DFd F p p<.05 ges
## 2 WMSize 2 30 8.199485 0.001442758 * 0.04500834
##
## $`Mauchly's Test for Sphericity`
## Effect W p p<.05
## 2 WMSize 0.03659746 8.793399e-11 *
##
## $`Sphericity Corrections`
## Effect GGe p[GG] p[GG]<.05 HFe p[HF] p[HF]<.05
## 2 WMSize 0.5093199 0.01137853 * 0.5113321 0.01128146 *
ezANOVA(data = Baypar_WMSize %>% filter(Exp =='Exp4a'), dv= InP_curve, wid=NSub, within = .(WMSize))
## Warning: Converting "NSub" to factor for ANOVA.
## $ANOVA
## Effect DFn DFd F p p<.05 ges
## 2 WMSize 2 30 1.507457 0.237775 0.01146743
##
## $`Mauchly's Test for Sphericity`
## Effect W p p<.05
## 2 WMSize 0.09784573 8.586046e-08 *
##
## $`Sphericity Corrections`
## Effect GGe p[GG] p[GG]<.05 HFe p[HF] p[HF]<.05
## 2 WMSize 0.5257197 0.2392211 0.5313463 0.2393698
ezANOVA(data = Baypar_WMSize %>% filter(Exp =='Exp4b'), dv= InP_curve, wid=NSub, within = .(WMSize))
## Warning: Converting "NSub" to factor for ANOVA.
## $ANOVA
## Effect DFn DFd F p p<.05 ges
## 2 WMSize 2 30 5.970484 0.00656511 * 0.08780338
##
## $`Mauchly's Test for Sphericity`
## Effect W p p<.05
## 2 WMSize 0.1083146 1.749077e-07 *
##
## $`Sphericity Corrections`
## Effect GGe p[GG] p[GG]<.05 HFe p[HF] p[HF]<.05
## 2 WMSize 0.5286291 0.02522631 * 0.5349067 0.0247749 *
Export data for spss
### change the table struction for spss
Bayparlist_list <-split(Baypar_WMSize %>% select(c("NSub","Exp","WMSize", "wp","InP_curve")), Baypar_WMSize$WMSize)
Bayparlist_spss = NULL
for (i in 1: length(Bayparlist_list)){
temp = Bayparlist_list[[i]]
WMSize = unique(temp$WMSize)
temp$WMSize = NULL
colnames(temp) = c('NSub','Exp', paste0('wp_',WMSize), paste0('InP_curve_',WMSize))
if(i == 1)
Bayparlist_spss = temp
else
Bayparlist_spss = left_join(Bayparlist_spss, temp, by=c("Exp", "NSub"))
}
write_csv(Bayparlist_spss, paste0(modelPath, '/rlt/Bayparlist_spss.csv'))
anova analysis
Anova on mean reproduction biases
mpredY_sub$RP_Bias = mpredY_sub$m_repDur-mpredY_sub$curDur
RP_bias_Anova <- ezANOVA(data = mpredY_sub, dv= RP_Bias, wid=NSub, within= .(curDur, WMSize), between = .(Exp), detailed = TRUE, return_aov = FALSE )
## Warning: Converting "NSub" to factor for ANOVA.
## Warning: "curDur" will be treated as numeric.
## Warning: Converting "Exp" to factor for ANOVA.
## Warning: The column supplied as the wid variable contains non-unique values
## across levels of the supplied between-Ss variables. Automatically fixing this by
## generating unique wid labels.
## Warning: There is at least one numeric within variable, therefore aov() will be
## used for computation and no assumption checks will be obtained.
RP_bias_Anova
## $ANOVA
## Effect DFn DFd SSn SSd F p
## 1 Exp 4 75 0.20632235 6.7357598 0.5743293 6.820927e-01
## 2 curDur 1 75 47.70152785 7.0760974 505.5914872 4.589992e-35
## 4 WMSize 2 150 0.04568339 0.6961599 4.9216486 8.506722e-03
## 3 Exp:curDur 4 75 0.94184534 7.0760974 2.4956695 4.985449e-02
## 5 Exp:WMSize 8 150 0.25075426 0.6961599 6.7536818 1.507627e-07
## 6 curDur:WMSize 2 150 0.11975287 0.3666986 24.4927739 6.240642e-10
## 7 Exp:curDur:WMSize 8 150 0.12549340 0.3666986 6.4167178 3.670623e-07
## p<.05 ges
## 1 0.013680912
## 2 * 0.762294526
## 4 * 0.003061808
## 3 * 0.059548049
## 5 * 0.016578279
## 6 * 0.007986470
## 7 * 0.008366110
# main effect of Duration F(1.177, 3.532) = 377.965, p < .001, ηp² = .863.
(RP_bias_Anova$ANOVA)$DFn[3] *(RP_bias_Anova$`Sphericity Corrections`)$GGe[1]
## numeric(0)
(RP_bias_Anova$ANOVA)$DFd[3] *(RP_bias_Anova$`Sphericity Corrections`)$GGe[1]
## numeric(0)
#Duration × Experiment, F(12, 240) = 2.506, p = .004, ηp² = .111
(RP_bias_Anova$ANOVA)$DFn[5] *(RP_bias_Anova$`Sphericity Corrections`)$GGe[2]
## numeric(0)
(RP_bias_Anova$ANOVA)$DFd[5] *(RP_bias_Anova$`Sphericity Corrections`)$GGe[2]
## numeric(0)
mpredY_sub <- ungroup(mpredY_sub)
res.aov <- rstatix::anova_test(data = mpredY_sub, dv = RP_Bias, wid = NSub, within = c(curDur, WMSize), between = Exp)
## Warning: The 'wid' column contains duplicate ids across between-subjects
## variables. Automatic unique id will be created
get_anova_table(res.aov, correction = "GG")
variance of prior (anova)
ezANOVA(data = Bayparlist, dv= sig_pr2_log, wid=NSub, between = .(Exp))
## Warning: Converting "NSub" to factor for ANOVA.
## Warning: Converting "Exp" to factor for ANOVA.
## Warning: The column supplied as the wid variable contains non-unique values
## across levels of the supplied between-Ss variables. Automatically fixing this by
## generating unique wid labels.
## Coefficient covariances computed by hccm()
## $ANOVA
## Effect DFn DFd F p p<.05 ges
## 1 Exp 4 75 0.5384571 0.70791 0.02791603
##
## $`Levene's Test for Homogeneity of Variance`
## DFn DFd SSn SSd F p p<.05
## 1 4 75 0.04608107 1.121874 0.7701578 0.5480179
variance of motor noise (anova)
ezANOVA(data = Bayparlist, dv= sig_mn2, wid=NSub, between = .(Exp))
## Warning: Converting "NSub" to factor for ANOVA.
## Warning: Converting "Exp" to factor for ANOVA.
## Warning: The column supplied as the wid variable contains non-unique values
## across levels of the supplied between-Ss variables. Automatically fixing this by
## generating unique wid labels.
## Coefficient covariances computed by hccm()
## $ANOVA
## Effect DFn DFd F p p<.05 ges
## 1 Exp 4 75 0.602902 0.6617231 0.03115306
##
## $`Levene's Test for Homogeneity of Variance`
## DFn DFd SSn SSd F p p<.05
## 1 4 75 0.001201703 0.03964056 0.568406 0.6863391
Ks anova
ezANOVA(data = Bayparlist, dv= ks, wid=NSub, between = .(Exp))
## Warning: Converting "NSub" to factor for ANOVA.
## Warning: Converting "Exp" to factor for ANOVA.
## Warning: The column supplied as the wid variable contains non-unique values
## across levels of the supplied between-Ss variables. Automatically fixing this by
## generating unique wid labels.
## Coefficient covariances computed by hccm()
## $ANOVA
## Effect DFn DFd F p p<.05 ges
## 1 Exp 4 75 85.93971 3.114241e-27 * 0.8208993
##
## $`Levene's Test for Homogeneity of Variance`
## DFn DFd SSn SSd F p p<.05
## 1 4 75 0.1584934 0.3319383 8.952723 5.765718e-06 *
mean of prior (anova)
ezANOVA(data = Bayparlist, dv= mu_pr, wid=NSub, between = .(Exp))
## Warning: Converting "NSub" to factor for ANOVA.
## Warning: Converting "Exp" to factor for ANOVA.
## Warning: The column supplied as the wid variable contains non-unique values
## across levels of the supplied between-Ss variables. Automatically fixing this by
## generating unique wid labels.
## Coefficient covariances computed by hccm()
## $ANOVA
## Effect DFn DFd F p p<.05 ges
## 1 Exp 4 75 0.9353417 0.4482419 0.04751463
##
## $`Levene's Test for Homogeneity of Variance`
## DFn DFd SSn SSd F p p<.05
## 1 4 75 0.02732829 1.139767 0.4495704 0.7723812
weight of prior
Weight of the prior \(w_p\)
Baypar_WMSize$wp = as.numeric(Baypar_WMSize$wp)
mwp <- Baypar_WMSize%>%dplyr::group_by(Exp, WMSize)%>% dplyr::summarise(m_wp = mean(wp), n= n(), se_wp = sd(wp)/sqrt(n-1))
## `summarise()` has grouped output by 'Exp'. You can override using the `.groups`
## argument.
plt_wp <- ggplot(mwp, aes(Exp, m_wp, ymin = m_wp - se_wp, ymax = m_wp + se_wp, group =interaction(Exp, WMSize), color = factor(WMSize)), shape = factor(WMSize))+
geom_line(stat = "identity",position = position_dodge(width = 0.2))+
geom_point(stat = "identity",position = position_dodge(width = 0.2))+
geom_errorbar(width=.2, position = position_dodge(width = 0.2)) +
#coord_cartesian(ylim = c(0.5, 1)) +
colorSet5+
labs(x = "", y = TeX("Weight of the prior $w_p$"), color = 'Memory Load') +
theme_new + theme(legend.position="top")
plt_wp
## geom_path: Each group consists of only one observation. Do you need to adjust
## the group aesthetic?
ggsave(paste0(getwd(), "/", modelPath, "/figures/plt_wp.png"), plt_wp, width = 3, height = 3)
## geom_path: Each group consists of only one observation. Do you need to adjust
## the group aesthetic?
combine InP and wp
fig4<-ggarrange(plt_wp, plt_InP, common.legend = TRUE, ncol=2, nrow=1, labels = c("a", "b"))
## geom_path: Each group consists of only one observation. Do you need to adjust
## the group aesthetic?
## geom_path: Each group consists of only one observation. Do you need to adjust
## the group aesthetic?
## geom_path: Each group consists of only one observation. Do you need to adjust
## the group aesthetic?
ggsave(paste0(getwd(), "/", modelPath, "/figures/fig4.png"), fig4, width = 6, height = 3)
fig4
fig5<-ggarrange(RP_bias_obs, plt_InP_linear, common.legend = TRUE, ncol=2, nrow=1, widths = c(5,3), labels = c("a", "b"))
## geom_path: Each group consists of only one observation. Do you need to adjust
## the group aesthetic?
ggsave(paste0(getwd(), "/", modelPath, "/figures/fig5.png"), fig5, width = 6, height = 3)
fig5
ANOVA analysis on wp
ezANOVA(data = Baypar_WMSize, dv= wp, wid=NSub, within=.(WMSize), between = .(Exp))
## Warning: Converting "NSub" to factor for ANOVA.
## Warning: Converting "Exp" to factor for ANOVA.
## Warning: The column supplied as the wid variable contains non-unique values
## across levels of the supplied between-Ss variables. Automatically fixing this by
## generating unique wid labels.
## $ANOVA
## Effect DFn DFd F p p<.05 ges
## 2 Exp 4 75 6.864458 9.177058e-05 * 0.26663773
## 3 WMSize 2 150 189.342145 9.270731e-42 * 0.01709301
## 4 Exp:WMSize 8 150 33.257671 1.198007e-29 * 0.01207081
##
## $`Mauchly's Test for Sphericity`
## Effect W p p<.05
## 3 WMSize 0.02558321 1.242606e-59 *
## 4 Exp:WMSize 0.02558321 1.242606e-59 *
##
## $`Sphericity Corrections`
## Effect GGe p[GG] p[GG]<.05 HFe p[HF] p[HF]<.05
## 3 WMSize 0.5064787 1.832257e-22 * 0.5067425 1.789171e-22 *
## 4 Exp:WMSize 0.5064787 4.021198e-16 * 0.5067425 3.954546e-16 *
ezANOVA(data = Baypar_WMSize%>%filter(Exp =='Exp2'), dv= wp, wid=NSub, within= .(WMSize) )
## Warning: Converting "NSub" to factor for ANOVA.
## $ANOVA
## Effect DFn DFd F p p<.05 ges
## 2 WMSize 2 30 73.10886 2.92466e-12 * 0.03584474
##
## $`Mauchly's Test for Sphericity`
## Effect W p p<.05
## 2 WMSize 0.03336632 4.604238e-11 *
##
## $`Sphericity Corrections`
## Effect GGe p[GG] p[GG]<.05 HFe p[HF] p[HF]<.05
## 2 WMSize 0.5084831 3.075084e-07 * 0.5103134 2.944652e-07 *
ezANOVA(data = Baypar_WMSize%>%filter(Exp =='Exp4a'), dv= wp, wid=NSub, within= .(WMSize) )
## Warning: Converting "NSub" to factor for ANOVA.
## $ANOVA
## Effect DFn DFd F p p<.05 ges
## 2 WMSize 2 30 45.80731 7.62135e-10 * 0.02711716
##
## $`Mauchly's Test for Sphericity`
## Effect W p p<.05
## 2 WMSize 0.01401012 1.059482e-13 *
##
## $`Sphericity Corrections`
## Effect GGe p[GG] p[GG]<.05 HFe p[HF] p[HF]<.05
## 2 WMSize 0.5035272 5.933032e-06 * 0.5042852 5.85177e-06 *
ezANOVA(data = Baypar_WMSize%>%filter(Exp =='Exp4b'), dv= wp, wid=NSub, within= .(WMSize) )
## Warning: Converting "NSub" to factor for ANOVA.
## $ANOVA
## Effect DFn DFd F p p<.05 ges
## 2 WMSize 2 30 72.16408 3.438001e-12 * 0.0546141
##
## $`Mauchly's Test for Sphericity`
## Effect W p p<.05
## 2 WMSize 0.02522245 6.493984e-12 *
##
## $`Sphericity Corrections`
## Effect GGe p[GG] p[GG]<.05 HFe p[HF] p[HF]<.05
## 2 WMSize 0.5063862 3.510759e-07 * 0.5077617 3.398985e-07 *
### pairwise.t.test on wp
pairwise.t.test(Baypar_WMSize$wp, Baypar_WMSize$Exp)
##
## Pairwise comparisons using t tests with pooled SD
##
## data: Baypar_WMSize$wp and Baypar_WMSize$Exp
##
## Exp1 Exp2 Exp3 Exp4a
## Exp2 2.5e-13 - - -
## Exp3 0.23 3.1e-08 - -
## Exp4a 0.23 4.1e-08 0.95 -
## Exp4b 0.95 2.0e-11 0.69 0.69
##
## P value adjustment method: holm
##independent T test
t.test((Baypar_WMSize%>%filter(Exp %in%c('Exp1')))$wp)
##
## One Sample t-test
##
## data: (Baypar_WMSize %>% filter(Exp %in% c("Exp1")))$wp
## t = 12.123, df = 47, p-value = 4.502e-16
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
## 0.2292634 0.3204897
## sample estimates:
## mean of x
## 0.2748766
t.test((Baypar_WMSize%>%filter(Exp %in%c('Exp2')))$wp)
##
## One Sample t-test
##
## data: (Baypar_WMSize %>% filter(Exp %in% c("Exp2")))$wp
## t = 15.978, df = 47, p-value < 2.2e-16
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
## 0.5204341 0.6703593
## sample estimates:
## mean of x
## 0.5953967
t.test((Baypar_WMSize%>%filter(Exp %in%c('Exp3')))$wp)
##
## One Sample t-test
##
## data: (Baypar_WMSize %>% filter(Exp %in% c("Exp3")))$wp
## t = 17.739, df = 47, p-value < 2.2e-16
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
## 0.3138573 0.3941508
## sample estimates:
## mean of x
## 0.3540041
t.test((Baypar_WMSize%>%filter(Exp %in%c('Exp4a')))$wp)
##
## One Sample t-test
##
## data: (Baypar_WMSize %>% filter(Exp %in% c("Exp4a")))$wp
## t = 11.932, df = 47, p-value = 7.948e-16
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
## 0.2968951 0.4173131
## sample estimates:
## mean of x
## 0.3571041
t.test((Baypar_WMSize%>%filter(Exp %in%c('Exp4b')))$wp)
##
## One Sample t-test
##
## data: (Baypar_WMSize %>% filter(Exp %in% c("Exp4b")))$wp
## t = 11.53, df = 47, p-value = 2.664e-15
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
## 0.2501705 0.3559257
## sample estimates:
## mean of x
## 0.3030481
standard variance of Ds \(\sigma_{s}\)
ezANOVA(data = Bayparlist, dv= sig_s2, wid=NSub, between = .(Exp))
## Warning: Converting "NSub" to factor for ANOVA.
## Warning: Converting "Exp" to factor for ANOVA.
## Warning: The column supplied as the wid variable contains non-unique values
## across levels of the supplied between-Ss variables. Automatically fixing this by
## generating unique wid labels.
## Coefficient covariances computed by hccm()
## $ANOVA
## Effect DFn DFd F p p<.05 ges
## 1 Exp 4 75 3.409153 0.01287802 * 0.1538485
##
## $`Levene's Test for Homogeneity of Variance`
## DFn DFd SSn SSd F p p<.05
## 1 4 75 0.01726871 0.1448521 2.235302 0.07314935
### pairwise.t.test on standard variance of $D_s$
pairwise.t.test(Bayparlist$sig_s2, Bayparlist$Exp)
##
## Pairwise comparisons using t tests with pooled SD
##
## data: Bayparlist$sig_s2 and Bayparlist$Exp
##
## Exp1 Exp2 Exp3 Exp4a
## Exp2 0.078 - - -
## Exp3 1.000 0.139 - -
## Exp4a 1.000 0.047 1.000 -
## Exp4b 1.000 0.013 1.000 1.000
##
## P value adjustment method: holm
##independent T test
t.test((Bayparlist%>%filter(Exp %in%c('Exp1')))$sig_s2)
##
## One Sample t-test
##
## data: (Bayparlist %>% filter(Exp %in% c("Exp1")))$sig_s2
## t = 7.2054, df = 15, p-value = 3.046e-06
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
## 0.02398402 0.04413435
## sample estimates:
## mean of x
## 0.03405919
t.test((Bayparlist%>%filter(Exp %in%c('Exp2')))$sig_s2)
##
## One Sample t-test
##
## data: (Bayparlist %>% filter(Exp %in% c("Exp2")))$sig_s2
## t = 3.0831, df = 15, p-value = 0.007574
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
## 0.0247191 0.1354437
## sample estimates:
## mean of x
## 0.08008139
t.test((Bayparlist%>%filter(Exp %in%c('Exp3')))$sig_s2)
##
## One Sample t-test
##
## data: (Bayparlist %>% filter(Exp %in% c("Exp3")))$sig_s2
## t = 11.682, df = 15, p-value = 6.232e-09
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
## 0.03170408 0.04585477
## sample estimates:
## mean of x
## 0.03877943
t.test((Bayparlist%>%filter(Exp %in%c('Exp4a')))$sig_s2)
##
## One Sample t-test
##
## data: (Bayparlist %>% filter(Exp %in% c("Exp4a")))$sig_s2
## t = 6.5979, df = 15, p-value = 8.467e-06
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
## 0.02042601 0.03992116
## sample estimates:
## mean of x
## 0.03017359
t.test((Bayparlist%>%filter(Exp %in%c('Exp4b')))$sig_s2)
##
## One Sample t-test
##
## data: (Bayparlist %>% filter(Exp %in% c("Exp4b")))$sig_s2
## t = 4.5291, df = 15, p-value = 0.0003995
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
## 0.01172269 0.03256513
## sample estimates:
## mean of x
## 0.02214391
# calculate conhens d
cohensD(sig_s2 ~ Exp,
data = Bayparlist%>%filter(Exp %in% c('Exp1', 'Exp2')),
method = "paired")
## Warning in cohensD(sig_s2 ~ Exp, data = Bayparlist %>% filter(Exp %in%
## c("Exp1", : calculating paired samples Cohen's d using formula input. Results
## will be incorrect if cases do not appear in the same order for both levels of
## the grouping factor
## [1] 0.4080594
cohensD(sig_s2 ~ Exp,
data = Bayparlist%>%filter(Exp %in% c('Exp2', 'Exp3')),
method = "paired")
## Warning in cohensD(sig_s2 ~ Exp, data = Bayparlist %>% filter(Exp %in%
## c("Exp2", : calculating paired samples Cohen's d using formula input. Results
## will be incorrect if cases do not appear in the same order for both levels of
## the grouping factor
## [1] 0.3966069
cohensD(sig_s2 ~ Exp,
data = Bayparlist%>%filter(Exp %in% c('Exp3', 'Exp4a')),
method = "paired")
## Warning in cohensD(sig_s2 ~ Exp, data = Bayparlist %>% filter(Exp %in%
## c("Exp3", : calculating paired samples Cohen's d using formula input. Results
## will be incorrect if cases do not appear in the same order for both levels of
## the grouping factor
## [1] 0.3980014
cohensD(sig_s2 ~ Exp,
data = Bayparlist%>%filter(Exp %in% c('Exp2', 'Exp4a')),
method = "paired")
## Warning in cohensD(sig_s2 ~ Exp, data = Bayparlist %>% filter(Exp %in%
## c("Exp2", : calculating paired samples Cohen's d using formula input. Results
## will be incorrect if cases do not appear in the same order for both levels of
## the grouping factor
## [1] 0.4654873
cv slope
ezANOVA(data = CVslopes, dv= obs_cv_slope, wid=NSub, between = .(Exp))
## Warning: Converting "NSub" to factor for ANOVA.
## Warning: Converting "Exp" to factor for ANOVA.
## Warning: The column supplied as the wid variable contains non-unique values
## across levels of the supplied between-Ss variables. Automatically fixing this by
## generating unique wid labels.
## Warning: Collapsing data to cell means. *IF* the requested effects are a subset
## of the full design, you must use the "within_full" argument, else results may be
## inaccurate.
## Coefficient covariances computed by hccm()
## $ANOVA
## Effect DFn DFd F p p<.05 ges
## 1 Exp 4 75 0.1126364 0.9777182 0.005971402
##
## $`Levene's Test for Homogeneity of Variance`
## DFn DFd SSn SSd F p p<.05
## 1 4 75 0.08653032 1.378199 1.17722 0.3277357
### pairwise.t.test on standard variance of $D_s$
pairwise.t.test(CVslopes$obs_cv_slope, CVslopes$Exp)
##
## Pairwise comparisons using t tests with pooled SD
##
## data: CVslopes$obs_cv_slope and CVslopes$Exp
##
## Exp1 Exp2 Exp3 Exp4a
## Exp2 1 - - -
## Exp3 1 1 - -
## Exp4a 1 1 1 -
## Exp4b 1 1 1 1
##
## P value adjustment method: holm
##independent T test
t.test((CVslopes%>%filter(Exp %in%c('Exp1')))$obs_cv_slope)
##
## One Sample t-test
##
## data: (CVslopes %>% filter(Exp %in% c("Exp1")))$obs_cv_slope
## t = -6.5917, df = 47, p-value = 3.406e-08
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
## -0.2970793 -0.1581479
## sample estimates:
## mean of x
## -0.2276136
t.test((CVslopes%>%filter(Exp %in%c('Exp2')))$obs_cv_slope)
##
## One Sample t-test
##
## data: (CVslopes %>% filter(Exp %in% c("Exp2")))$obs_cv_slope
## t = -7.47, df = 47, p-value = 1.59e-09
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
## -0.3369816 -0.1939864
## sample estimates:
## mean of x
## -0.265484
t.test((CVslopes%>%filter(Exp %in%c('Exp3')))$obs_cv_slope)
##
## One Sample t-test
##
## data: (CVslopes %>% filter(Exp %in% c("Exp3")))$obs_cv_slope
## t = -6.932, df = 47, p-value = 1.037e-08
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
## -0.3153826 -0.1735021
## sample estimates:
## mean of x
## -0.2444424
t.test((CVslopes%>%filter(Exp %in%c('Exp4a')))$obs_cv_slope)
##
## One Sample t-test
##
## data: (CVslopes %>% filter(Exp %in% c("Exp4a")))$obs_cv_slope
## t = -5.4143, df = 47, p-value = 2.047e-06
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
## -0.3547829 -0.1625578
## sample estimates:
## mean of x
## -0.2586703
t.test((CVslopes%>%filter(Exp %in%c('Exp4b')))$obs_cv_slope)
##
## One Sample t-test
##
## data: (CVslopes %>% filter(Exp %in% c("Exp4b")))$obs_cv_slope
## t = -7.4396, df = 47, p-value = 1.768e-09
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
## -0.3468257 -0.1991797
## sample estimates:
## mean of x
## -0.2730027
# calculate conhens d
cohensD(obs_cv_slope ~ Exp,
data = CVslopes%>%filter(Exp %in% c('Exp1', 'Exp2')),
method = "paired")
## Warning in cohensD(obs_cv_slope ~ Exp, data = CVslopes %>% filter(Exp %in% :
## calculating paired samples Cohen's d using formula input. Results will be
## incorrect if cases do not appear in the same order for both levels of the
## grouping factor
## [1] 0.1178123
cohensD(obs_cv_slope ~ Exp,
data = CVslopes%>%filter(Exp %in% c('Exp2', 'Exp3')),
method = "paired")
## Warning in cohensD(obs_cv_slope ~ Exp, data = CVslopes %>% filter(Exp %in% :
## calculating paired samples Cohen's d using formula input. Results will be
## incorrect if cases do not appear in the same order for both levels of the
## grouping factor
## [1] 0.06429713
cohensD(obs_cv_slope ~ Exp,
data = CVslopes%>%filter(Exp %in% c('Exp3', 'Exp4a')),
method = "paired")
## Warning in cohensD(obs_cv_slope ~ Exp, data = CVslopes %>% filter(Exp %in% :
## calculating paired samples Cohen's d using formula input. Results will be
## incorrect if cases do not appear in the same order for both levels of the
## grouping factor
## [1] 0.04045837
cohensD(obs_cv_slope ~ Exp,
data = CVslopes%>%filter(Exp %in% c('Exp2', 'Exp4a')),
method = "paired")
## Warning in cohensD(obs_cv_slope ~ Exp, data = CVslopes %>% filter(Exp %in% :
## calculating paired samples Cohen's d using formula input. Results will be
## incorrect if cases do not appear in the same order for both levels of the
## grouping factor
## [1] 0.01404826
Model prediction error
m_predErr_sub<- mpredY_sub%>%
dplyr::group_by(Exp, WMSize, NSub) %>% dplyr::summarise(
mpredRP_err=mean(predRP_err),
mpredVar_err=mean(predVar_err),
mpredcv_err = mean(predcv_err),
mpredRP_rerr = mean(predRP_rerr),
mpredVar_rerr = mean(predVar_rerr),
mpredcv_rerr = mean(predcv_rerr))
## `summarise()` has grouped output by 'Exp', 'WMSize'. You can override using the
## `.groups` argument.
m_predErr<- m_predY%>%
dplyr::group_by(Exp, WMSize) %>% dplyr::summarise(
mmpredcv_err = mean(mpredcv_err),
mmpredRP_err=mean(mpredRP_err),
mmpredVar_err=mean(mpredVar_err),
mmpredRP_rerr = mean(mpredRP_rerr),
mmpredVar_rerr = mean(mpredVar_rerr),
mmpredcv_rerr = mean(mpredcv_rerr))
## `summarise()` has grouped output by 'Exp'. You can override using the `.groups`
## argument.
m_predErr
plt_rErrorScatter = ggplot(m_predErr_sub, aes(mpredRP_rerr*100, mpredcv_rerr*100, color = WMSize, alpha = .9)) +
#geom_hline(yintercept = 0, linetype='dashed')+ geom_vline(xintercept = 0, linetype='dashed')+
geom_point() +
geom_point(data = m_predErr, aes(mmpredRP_rerr*100, mmpredcv_rerr*100, color = WMSize, alpha = .9, size = 1 ))+
xlab('Relative prediction error of reproduction mean(%)')+ ylab('Relative prediction error of CV (%)')+colorSet3+
facet_wrap(~Exp)+
theme_new+ theme(legend.position = 'top')+guides(size="none")+guides(alpha="none")
plt_rErrorScatter
plt_rErrorScatter = ggplot(m_predErr_sub, aes(mpredRP_rerr*100, mpredVar_rerr*100, color = WMSize, alpha = .9)) +
#geom_hline(yintercept = 0, linetype='dashed')+ geom_vline(xintercept = 0, linetype='dashed')+
geom_point() +
geom_point(data = m_predErr, aes(mmpredRP_rerr*100, mmpredVar_rerr*100, color = WMSize, alpha = .9, size = 1 ))+
xlab('Relative prediction error in the RP means (%)')+ ylab('Relative prediction error in the RP variance (%)')+colorSet3+
facet_wrap(~Exp)+
theme_new+ theme(legend.position = 'top')+guides(size="none")+guides(alpha="none")
plt_rErrorScatter
plt_ErrorScatter = ggplot(m_predErr_sub, aes(mpredRP_err, mpredVar_err, color = WMSize, alpha = .9)) +
geom_hline(yintercept = 0, linetype='dashed')+ geom_vline(xintercept = 0, linetype='dashed')+
geom_point() +
geom_point(data = m_predErr, aes(mmpredRP_err, mmpredVar_err, color = WMSize, alpha = .9, size = 1 ))+
xlab('Prediction error in the RP means (s)')+ ylab('Prediction error in the RP variance (s)')+colorSet3+
facet_wrap(~Exp)+
theme_new+ theme(legend.position = 'top')+guides(size="none")+guides(alpha="none")
ggsave(paste0(getwd(), "/", modelPath, "/figures/plt_ErrorScatter.png"), plt_ErrorScatter, width = 7, height = 4)
plt_ErrorScatter
model comparison (logarithmic vs. linear)
m_predErr_sub$model = 'logarithmic'
m_predErr$model = 'logarithmic'
linear_model = 'linear_rstan'
m_predErr_linear = read.csv(paste0(getwd(), "/", rstanmodelPath, '/models/', linear_model, "/rlt/m_predErr_", linear_model, ".csv"))
m_predErr_linear$X = NULL
m_predErr_sub_linear = read.csv(paste0(getwd(), "/", rstanmodelPath, '/models/', linear_model, "/rlt/m_predErr_sub_", linear_model, ".csv"))
m_predErr_sub_linear$X = NULL
m_predErr_sub_all = rbind(m_predErr_sub, m_predErr_sub_linear)
m_predErr_all = rbind(m_predErr, m_predErr_linear)
m_predErr_all$WMSize = as.factor(m_predErr_all$WMSize)
levels(m_predErr_all$WMSize) = c("low", "medium", "high")
temp = m_predErr_all %>% filter(model == 'logarithmic') %>%summarise(abs_mmpredcv_err = abs(mmpredcv_err))
plt_Err_CV_all = ggplot(m_predErr_all, aes(abs(mmpredRP_err), abs(mmpredcv_err), group = interaction(model, Exp, WMSize), color = model, shape = Exp, size = WMSize)) +
geom_point() +
geom_hline(yintercept = round(max(temp$abs_mmpredcv_err), 4)+0.0005, linetype='dashed')+
xlab('Prediction error in the RP means (s)')+ ylab('Prediction error in CV')+colorSet3+
scale_shape_manual(values = c(6, 7, 16, 17,8)) +
theme_new+
theme(legend.position = 'top')+
labs(size = 'Memory Load')+
guides(colour = guide_legend(order = 1, nrow=2,byrow=TRUE),
shape = guide_legend(order =2, nrow=2,byrow=TRUE),
size = guide_legend(order = 3, nrow=3,byrow=TRUE))
ggsave(paste0(getwd(), "/", modelPath, "/figures/plt_Err_CV_all.png"), plt_Err_CV_all, width = 7, height = 4)
## Warning: Using size for a discrete variable is not advised.
plt_Err_CV_all
## Warning: Using size for a discrete variable is not advised.
m_predY_acc = m_predErr_sub_all%>%
dplyr::group_by(Exp, model) %>%
dplyr::summarize(mmpredRP_rerr = mean(mpredRP_rerr)*100,
mmpredVar_rerr = mean(mpredVar_rerr)*100,
mmpredcv_rerr = mean(mpredcv_rerr)*100,
mmpredRP_acc = (1-mean(mpredRP_rerr))*100,
mmpredVar_acc = (1-mean(mpredVar_rerr))*100,
mmpredCV_acc = (1-mean(mpredcv_rerr))*100)
## `summarise()` has grouped output by 'Exp'. You can override using the `.groups`
## argument.
m_predY_acc